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Abstract 

For  precise  time  intercomparisons  between  a  master  frequency  standard  and  a  slave  time  scale,  we 
fiavc/ountl  it  useful  to  quantitatively  compare  different  fittinfi  strategies  by  examining  the  standard 
uncertainty  in  time  or  average /requeue}?.  It  is  particularly  useful  when  designing  procedures  which 
use  intermittent  intercomparisons,  with  some  parameterized  fit  used  to  interpolate  or  extrapolate 
from  the  calibrating  intercomparisons.  He  use  the  term  '*metafi(ting"for  the  choices  that  are  made 
before  a  fitting  procedure  is  operationally  adopted.  He  present  methods  for  calculating  the  standard 
uncertainty  for  general,  weighted  least-squares  fits  and  a  method  for  optimizing  these  weights  for  a 
general  noise  model  suitable  for  many  PTTI  applications.  He  present  the  results  of  the  metafitting 
of  procedures  for  the  use  of  a  regular  schedule  of  (hypothetical)  high-accuracy  frequency  calibration 
of  a  maser  time  scale.  We  hove  identified  a  cumulative  series  of  improvements  that  give  a  significant 
reduction  of  the  expected  standard  uncertainty,  compared  to  the  simplest  procedure  of  resetting 
the  maser  synthesiser  after  each  calibration.  The  metafitting  improvements  presented  include  the 
optimum  choice  of  weights  for  the  calibration  runs,  optimised  over  a  period  of  a  week  or  10  days. 


Introduction 

(n  preparing  to  fit  precision  time  compari.son  data,  usually  questions  concerning  “optimal”  fitting 
.strategies  have  been  addres.sed  in  a  generic  rather  than  in  a  specific  sense.  It  is  interesting  to 
examine  whether,  for  specific  cases,  significant  advantages  might  accrue  from  customizing  the  fitting 
strategy  to  the  specific  pattern  of  data  points  and  the  noise  spectrum.  In  practice,  many  really 
important  choices  are  made  before  any  fit  is  finalized,  and  yet  are  not  necessairly  optimized  as  part 
of  the  fitting  procedure,  (a)  A  fitting  metric  and  method  must  be  chosen  (such  as  Jeast-,squares 
fitting).  (1))  The  set  of  parameterized  basis  function, s  mu.st  be  chasen;  basis  function  number  and 
type  (such  as  a  second  order  polynomial),  (c)  An  outlier  removal  method  may  be  adopted  (such 
as  iteratively  discarding  a  limited  number  of  points  having  anomalously  high  residuals),  (d)  The 
relative  weighting  to  be  given  to  each  data  point  must  be  determined  (such  as  the  use  of  end-point 
only  linear  fits  vs  unweighted  linear  least-squares  fits),  (e)  A  final  “consistency  of  fit  with  data 
and  noise  model”  parameter  should  be  derived  (such  as  the  reduced  for  a  least-squares  fit  with 
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known  noise).  Many  of  these  choice.s  depend  in  subtle  ways  on  the  type  of  noise  encountered,  and 
here  precise  time  measurements  often  provide  details  of  the  noise  .spectrum  which  are  not  trivial  to 
incorporate  into  an  optimal  treatment  of  (a)-(e).  Since  the  term  “fitting”  is  generally  interpreted 
as  referring  to  the  determination  of  a  set  of  parameters  from  a  particular  data  set,  we  use  the  term 
“metafitting”  to  encompass  the  optimization  of  the  broader  processes  such  as  (a)-(e). 

But  in  what  sense  is  this  “metafitting”  to  be  judged?  At  first  sight,  there  appear  to  be  too  many 
choices.  The  fitting  Bright  be  optimized  in  an  average  .sense,  minimizing  some  metafitting  metric, 
function  that  sums  over  experimental  residuals.  If  the  autocorrelation  function  of  the  noise  is 
known  (or  modeled)  it  is  possible  to  calculate  and  minimize  the  metric  function  summing  over  the 
expected  “residuals”  at  unmeasured  times.  Thus  the  fitting  might  be  optimized  in  a  local  sense, 
minimizing  a  residual  at  a  specific  time  (open  to  choice),  or  it  might  be  optimized  to  minimize 
a  residual  of  the  average  frequency  over  a  specific  interval  (each  end  being  open  to  choice).  The 
point,  points  or  interval  must  be  cho.sen,  and  a  procedure  must  be  found  to  estimate  the  expected 
residual(s)  at  times  other  than  those  at  which  measurements  have  Ijeen  taken. 

Fortunately,  international  guidelines  [1]  now  strongly  .suggest  a  good  quantity  to  optimize:  the 
“standard  uncertainty”,  which  is  the  root-mean-square  residual  of  the  fit’s  extrapolation  or  inter¬ 
polation  to  a  specific  point,  one  not  necassairly  included  in  the  fit.  It  is  al.so  a  good  quantity 
to  optimize  in  that  a  standard  method  [2]  (the  Wiener- Kolmogoroff  theory)  exists  for  any  fitting 
procedure  that  uses  a  linear  combination  of  the  data.  All  least-squares  fitting  procedures  with 
linear  coefficients  can  be  handled  explicitly  in  this  way  [5|.  For  the  purposes  of  time  and  frequency 
metrology,  metafitting  to  minimize  the  standard  uncertainty  is  a  good  choice  -  but  it  might  not 
be  as  good  a  choice  in  other  applications  (where,  for  example  it  might  be  more  appropriate  to  try 
minimizing  the  occurrence  of  outlier  events  having  disasterous  consequences).  For  frequency  or 
time  interval  metrology,  the  standard  uncertainty  in  the  average  frequency  over  an  interval  makes 
an  even  more  attractive  discriminant  for  metafitting. 

The  pow'er  law  noise  models  appropriate  for  PTTI  phase  comparisons  can  have  low  frequency 
divergences  that  appear  to  be  worrying  to  some  purists  who  wish  to  assure  strict  stationarity  of 
any  proce.ss  before  developing  its  formalism.  In  the  development  of  computable  forms  [3], [5],  it  is 
straightforward  to  show  that  the  standard  uncertainty  in  the  average  frequency  of  a  least-.squares 
fit  is  not  divergent  for  most  commonly  encountered  power  law  noise  spectra,  with  the  exception  of 
random  walk  frequency  noise.  However,  the  real  question  is  more  stringent  than  simple  stationarity; 
have  we  enough  long-term  data  on  the  system  being  modeled  to  obtain  results  which  converge? 
We  believe  that  this  type  of  question  can  Ije  rigorously  handled  by  imposing  a  low-frequency  cutoff 
(and  thus  ensuring  a  formal  stationarity),  and  then  verifying  not  merely  that  any  re.sults  extracted 
from  the  model  converge  as  the  low  frequency  limit  approaches  zero  -  but  also  that  the  results  have 
converged  to  the  desired  degree  before  the  low  frequency  limit  is  sampling  Fourier  components  of 
the  noise  which  have  not  been  measured. 


Choosing  Weights  in  Weighted  Least-Squares  Fits 


We  present  here  a  general  strategy  for  evaluating  and  optimizing  distriliutions  of  weights  in  a 
weighted  least-squares  fit  to  phase  data.  We  will  concentrate  on  optimizing  fitting  that  compares  the 
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frequency  of  a  continuously  operated  oscillator  with  a  frequency  standard  (perhaps  intermittently 
operated),  for  the  purposes  of  frequency  calibration.  The  strategy  is  based  upon  the  analytic 
expressions  for  the  standard  uncertainty  in  frequency,  generally  extrapolated  over  a  wider  interval 
than  the  calibration  interval,  where  a  dense  set  of  high-precision  phase  comparisons  would  normally 
be  available.  A  noise  model  is  assumed  which  has  wide  applicability  to  a  broad  range  of  frequency 
standards.  The  degree  of  frequency  control  will  be  evaluatalrie  for  any  set  of  weights  in  a  weighted 
least-squares  fit  that  is  linear  in  the  fitting  coefficients  (Irut  fully  general  in  the  choice  of  Irasis 
functions). 

As  our  metafitting  metric  we  choose  the  standard  uncertainty  in  the  average  frequency,  evaluated 
over  a  general  interval  which  could  be  considerably  displaced  from  the  fitting  interval.  This  is 
the  most  appropriate  metric  for  frequency  metrology  applications,  since  the  standard  uncertainty 
is  now  the  internationally  recommended  [1]  way  of  specifying  calibration  uncertainty.  With  our 
procedures  the  standard  uncertainty  in  average  frequency  can  be  evaluated  for  a  broad  class  of 
noise  models,  for  any  set  of  fitting  points,  for  any  extrapolation  or  interpolation  interval,  for  any 
linear  combination  of  arl)itrary  basis  functions,  and  for  any  least-squares  weighting.  In  particular, 
in  any  of  the  above  cases  our  procedure  can  evaluate  the  standard  uncertainty  the  equal-weight 
procedure  (advocated  for  its  robustness)  to  the  end-point  procedure  (advocated  for  its  “optimum” 
estimate  of  frequency  for  some  pure  c]as.ses  of  noise),  as  well  as  any  intermediate  case  with  higher 
weights  near  the  end-points  of  the  calibration  interval.  The  procedure  permits  the  evaluation  of  the 
trade-off  of  uncertainty  for  other  proc^edures  which  are  percieved  as  being  more  robust.  As  is  shown 
below,  even  with  large  data  sets,  in  some  cases  it  appears  to  be  feasible  to  choose  the  optimum 
set  of  weights  which  minimize  the  standard  uncertainty  in  average  frequency  for  the  interval  being 
considered. 


Noise  Model 


The  noise  model  im(f)  is  the  modeled  phase  difference  Iretween  the  master  frequency  standard 
and  the  standard  being  calibrated.  The  noise  mode]  is  taken  as  being  the  sum  of  a  deterministic 
part  (which  could  include  a  phase  offset,  frequency  offset  and  frequency  drift)  and  a  random  noise 
part,  a;o(t).  The  random  noise  includes  the  “full”  noise  model  that  is  usually  used  in  discussions  of 
frequency  standard  stability  [9]:  a  sum  of  five  noise  processes,  each  normally  distributed  about  the 
mean  (but  with  variances  which  depend  on  the  time  sampled  in  different  ways)  that  have  spectral 
densities  of  phase  noise  that  are  power  laws  which  range  from  flat  to  increasingly  divergent 

at  low  frequencies.  Exj>ressing  the  five  terms  in  terms  of  the  spectral  density  of  the  mean-square 
of  the  fluctuations  in  (or  yo(0)  a  frequency  /,  5’y(/),  each  noise  term  is  described  by  an 

amplitude  ha  which  is  taken  to  be  independent  of  any  time  translations  (station arity  and  random 
phase  approximations).  The  sum  includes  a  =  2,  white  phase  noise  in  x;  cv  =  1,  flicker  (1//)  noise 
in  x;  fv  =  0,  white  frequency  noise  and  random  walk  phase  noise;  a  =  —  1,  flicker  frequency  noise; 
and  a  —  —2,  random  walk  frequency  noise.  A  low-frequency  cutoff  fi  and  an  upper  frequency  cutoff 
fh-  The  spectral  density  of  the  mean-square  fluctuations  in  xo(t)  is  Sx{f),  and  for  this  noise  model 
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(I) 


s,if)  =  i;  Kr  and  SM)  =  E 

For  a  given  noise  model  of  this  type,  the  standard  uncertainty  cjf  the  fit  at  any  given  time  can  lie 
calculated  from  the  autocorrelation  function  (a;o(t)2:o(^  +  "J"))-  It  is  divergent  for  four  of  our  five 
types  of  noise  unless  a  low-frequency  cutoff  is  applied,  and  even  then  can  challenge  the  accuracy 
and  dynamic  range  capacities  of  classical  computing.  Analytic  expressions  for  this  autocorrelation 
function  exist  for  each  type  of  noise  [5],  and  modern  arbitrary-precision  computer  languages  are 
able  to  cope  directly  with  the  autocorrelation  function. 

In  our  analysis  of  the  uncertainty  associated  with  any  useful  lea.st-squares  fit,  we  expect  no  di¬ 
vergences  to  infinity  in  the  standard  uncertainty,  and  so  the  combinations  of  the  autocorrelation 
functions  must  have  their  divergent  parts  cancel,  with  the  fitting  itself  acting  as  low-frequency  cut¬ 
off.  In  considering  the  standard  uncertainty  of  average  frequency  from  a  least-squares  fit,  we  have 
found  it  helpful  to  use  analytic  expressions  [5],  [4j,  [If]  for  the  less  divergent  general  two-interval 
covariance  of  the  random  noise  model,  that  is  the  covariance  of  the  time-scale  departure  over  the 
time  interval  [ti,t2]  with  the  time-scale  departure  over  the  time  interval 


s  =  ([xo(f.2)  -  XQ{ti)]  [0:0(4)  -  ;co(f3)]) 

=  r  f  \yo{nyo(n)dt"dt'  (2) 

■■'q  Jh 

=  {(-2  ~  h){h  -  h)  (y[qy2iy[t:,,t4l) 


where  is  the  general  covariance  of  the  average  frequency:  a  generalization  of  the 

twv>-.sample  variance  of  the  average  frequency.  The  generalization  includes  the  possibility  of  an 
overlap  of  the  intervals  (as  w'ell  as  the  possibility  of  a  "dead  time”  between  the  intervals),  and 
incorporates  the  possil)ility  of  considering  the  frequency  average  over  tw'o  time  intervals  of  different 
duration.  Just  as  for  the  two-sample  variance  of  y,  and  for  the  autocorrelation  function  of  x(t),  the 
covariance  separates  into  the  five  terms  of  the  noise  model. 

Analytic  fornus  for  the  five  terms  of  the  autocorrelation  function  of  x{t)  and  for  the  five  terms 
of  the  general  cross-correlation  of  y  are  given  in  references  [5],  [4]  and  [3],  derived  with  only  the 
usual  assumptions  about  high  and  low  frequency  limits  to  the  noise  bandwidth.  The  references  also 
contain  some  comments  on  prachical  methods  for  computing  values  using  these  forms. 


Weighted  Least-Squares  Fits 


Weighted  least-squares  fitting  cJiooses  the  n  linear  coeffic.ients  dj  of  the  n  basis  functions  gi(t), 
to  arrive  at  a  function  Xp{t.)  w'hich  w’ill  be  used  for  interpolation  or  extrapolation.  In  frequency 
standards  work,  w'e  would  usually  fit  a  phase  offset,  a  frequency  offset,  and  sometimes  a  drift  rate 
and  higher  terms  such  as  daily  or  seasonal  fluctuations. 
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(3) 


Xp{t)  =  d  ■  g{t.)  =  dj  +  d^t  +  d^l^  +  ^  digi{<). 

1=4 


The  coefficient  vector  (Tis  chosen  to  minimize  the  sum  over  the  A'  fitting  points  with  phase  difference 
values  of  x(t.i)  at  times 


~  9{u.)f  (4) 

i=l 

where  the  weight  VKj  is  applied  to  the  square  of  the  residual.  Least-squares  fitting  is  done  by 
setting  the  n  derivatives  of  C2  to  zero,  which  gives  a  set  of  n  linear  equations  which  can 

be  solved  for  the  n  fitting  coefficients  of  d:  Gd  =  where  G  is  an  n  x  n  matrix  with  elements 
Cfjr  —  ^yifjq{k)Sr{k),  and  s  is  an  n-dirnensional  vectoz'  with  elements  Sr  — 

For  the  purposes  of  modelling  the  standard  uncertainty,  we  usex(j(k)  to  model  since  it  can  be 
shown  [5]  that  any  general  offset  in  phase,  offset  in  frefjuency  or  a  linear  frequency  drift  is  exactly 
absorbed  by  the  fit. 


Metafitting  with  Time  Uncertainty  Metric 


One  candidate  metric  for  judging  weighted  least-squares  fits  is  the  standard  uncertainty  in  time, 
determined  at  a  specific  time  f,  relative  to  the  set  of  fitting  poiirts  {/,,).  \Vc  can  explicitly  calculate 
the  effects  of  the  weighted  least-squar’es  fit  reacting  to  tfie  noise  mode!  for  this  time  we  arc  not 
lestricted  to  studying  the  variance  at  the  fitting  points.  The  expected  variance  in  x(t)  from  the  fit 
d-  s{t)  can  be  calculated  in  terms  of  the  autocorrelation  function  Ic.ft  <  X(,{ti)x(](lj), 


x{t)  -  d-g{t) 


N  N 
i=0j=0 


(5) 


w’here  j[>o(0  =  1  and  hot  the  standard  noise  model, 

the  autocorrelation  function  {xo{t)xci{t))  can  be  evaluatecf  analytically  [5j,  although  the  result¬ 
ing  expressions  can  challenge  the  dynamic  range  of  conventional  computing.  The  square  root  of 
this  variance  in  2;(t)  would  be  the  formal  metric.  The  minimization  problem,  for  optimizing  this 
iiretric.  with  respect  to  the  weights  VF,-,  looks  intractable,  but  foi'  cases  of  most  interest  it  can  be 
substantially  simplified  in  the  same  way  as  is  described  below^  for  the  frequency  uncertainty  metric. 

Variants  of  this  L2  inelafitting  metric  are  also  possible,  summing  variances  over  multiple  test  times. 
Other  metafitting  metrics  of  the  Lp-norm  (Holder  norm)  class,  could  also  be  constructed.  The 
min- max  (limp  — ^  00)  norm  w'ould  ininiiiiize  the  maximum  expected  time  deviation  amongst  the 
test  times.  Metafitting  with  the  p  ~  1  metric  would  (for  this  class  of  metrics)  give  the  most  leeway 
in  allowing  a  small  number  of  test  points  to  have  large  variances.  All  these  metafitting  variants  arc 
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substantially  more  intricate  to  use,  and  do  not  readily  yield  the  major  computational  simplifications 
which  can  l>e  found  for  the  single- point  L2  metric. 

The  method  outlined  above  does  not  bring  any  great  new  insights  into  optimal  ways  of  combining 
equivalent  clocks,  nor  for  the  optimal  use  of  continuously  operated  primary  standards,  however  when 
a  .secondary  time  scale  is  to  track  a  primary  time  scale  where  only  intermittent  intercomparisons 
are  availaljle,  an  optimal  choice  could  be  made  in  terms  of  the  noise  processes  known  to  be  present. 


Metafitting  with  Frequency  Uncertainty  Metric 


For  precise  time  interval  work,  where  the  average  frequency  is  the  chief  quantity  of  intere.st,  we 
wish  to  minimize  the  standard  uncertainty  in  average  frequency  over  an  interval  [t,  f  -H  r],  caased 
by  the  noise  model  as  filtered  by  the  weighted  least-squares  fitting  procedure  to  the  points  {t,:}. 
Although  the  noise  model  is  independent  of  time  translations,  clearly  the  standard  uncertainty  in 
average  frequency,  would  be  expected  to  depend  on  the  offset  of  t.  from  {t,;},  as  well  as  the 
interval  breadth  r.  It  is  defined  by 


uJ(t,T) 


T 


} 


d-  {9{t  +  r)  -  g{t)] 


(6) 


We  note  that  {a:o(^+r)  —  2:Q(f)}  =  — 2;o(f7-i)]j  if  we  define  tj-o  =  t  and  —  t  +  r. 

Although  it  might  be  convenient  to  envisage  the  set  of  {tj)  as  an  ordered  set  with  tj  > 
it  is  not  necessary  to  do  so.  Ordering  the  fitting  points  does  not  detract  from  the  generality  in 
any  way,  but  we  do  not  wish  to  restrict  the  values  of  t.  or  t  -t-  r.  We  would  like  to  re-expre.ss  the 
d  ■  {git  +t)  —  as  a  sum  over  only  differences  of  the  form  xo(ti}  —  Xoitj).  We  note  that  we  can 
expand  xo(f.i)  =  -  xo(f,-i)},  so  that  d'  {ff(t  +  r)  -  is  equal  to 


G-'  ^  IVixo(fi)  ■  {g(t  +  r)  -  g(t)}  = 

^  m^l^oitj)  -  xo(tj-2)}  X^(G~^)^rgr(ti){gq(t  +  t)  -  gq(f)}  (7) 

i=l  j=2  g=Jr=l 

+2ro(fl)  53  53  ^(^~^)^rgr(t^){g^(t  +  r)  - 

8=1  1  r=] 

and  the  last  term,  multiplying  a:o(fi),  can  be  shown  to  ]>e  equal  to  zero.  To  show  this,  it  is  sufficient 
to  show  that  ^ig(ti)G~^g(t)  is  independent  of  f,,  or  that  is  equal  to  the 

vector  [1,0,0,...0].  We  oUserve  that,  from  our  definition  of  G  and  .since  gi  =  1,  G[1,0,0, ...0]  = 
and  premultiplying  by  G”'^  completes  this  proof,  provided  only  that  <71  is  a  constant. 
Thus  Uy(Tr)r^  is  equal  to 
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>  = 


(8) 


N+\ 

<  [X!  {^o(fi) 

.)=i 

-  X  X  X^^  ^)gr9r{h){9q{^  +  '’')“  5<j(f')}P 

i=l  j=2  <7=1  r=l 

M  n  n 

<  [{«o(iiV+i)  -  a;o(^Q)}  -  X  ~  xo(ti-i)}  EE(G'‘  }qr9fiti){9qi^  59(0}]  ^ 

i=l  q=lr=l 

Collecting  the  expressions  with  the  same  difference  term  {^o(t/)  -  xo(tj_j))  allows  us  to  write  a 
useful  form,  namely 


A^+l 

=<  [X  -  xo(tj_i))p  >,  (9) 

where  for2  <  j  <  iV,  Dj{t,T)  =  = 

1  and  Dj=/v>i(f.,r)  =  1.  Multiplying  the  terms  explicitly  gives  a  computable  form  for  the  standard 
uncertainty  in  average  frequency; 

y+iAf+i 

X  X  r)Dkit,T)  <  [xoitj)  ~  Xo(t7-l)][Xo(t;r)  -  Xo(tfc_l)]  >  .  (10) 

,7  =  1  fc=l 

The  utility  of  this  form  lies  in  the  fact  that  it  is  a  sum  over  functions  of  the  general  form  of  Eq.  2, 
which  ai'e  easier  to  compute  for  our  full  noise  model. 


Metafitting  Weights  for  Large  Data  Sets 


For  a  given  noise  model  (defined  by  the  5  parameters  {/laj  used  to  define  Sy(f),  and  a  given 
distribution  of  fitting  points  and  for  a  given  interval  [/, /.  +t];  the  standard  uncertainty  in 
average  frequency  over  the  interval  can  be  calculated:  u^(t,r).  Thus  a  choice  of  weights  can  be 
determined  which  minimizes  Uy(t.,T),  the  standard  uncertainty  due  to  the  effects  of  the  random 
noi.se.  For  each  fitting  point  added,  another  weight  must  be  determined.  For  small  sets  of  fitting 
points,  the  minimization  problem  is  tractable,  but  for  larger  sets  the  minimization  appears  much 
less  straightforward.  The  weights  could  be  parameterized  to  reduce  the  dimensionality  of  the 
problem,  at  the  expense  of  generality. 

The  full  generality  can  be  retained  by  largely  linearizing  the  problem.  For  N  fitting  points,  there 
are  also  N  weights  to  choose.  Without  loss  of  generality,  the  set  of  weights  {W;}  can  be  normalized: 

Wi  =  1.  If  the  partial  derivative  of  with  respec’t  to  can  be  constrained  to  be  zero,  then 
most  of  the  iV-dimensionaJ  search  problem  can  !>e  linearized,  lea\'ing  a  nonlinear  search  over  at  worst 
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[(n(n  +  l)/2)  -  1]  dimensions.  G~'  will  Ije  independent  of  Wk  if  eac;h  element  of  G  is  constrained 
to  be  a  constant,  Gqr  =  I2iLi  Since  G^r  ~  G,q,  and  since  normally  gi  =  1,  there 

remain  [(n.(n,  +  l)/2)  —  1]  values.  These  constraint  equations  are  used  in  the  linear  solution,  and 
the  optimum  values  of  Gqr  can  subsequently  be  found  by  nonlinear  searching  techniques. 

For  polynomial  fitting,  with  n  basis  functions  {gk(f}  =  the  partial  derivative  of  G“^  with 

respect  to  [N  —  2n  +  2]  VF,:’s  there  would  be  only  [2n  —  2]  dimensions  for  the  non-linear  search,  and 
if  the  problem  can  be  set  up  symmetrically  about  the  time  origin,  so  that  the  first  moment  of  the 
weights  and  all  odd  moments  are  zero,  there  would  be  only  [n  —  Ij  non-linear  search  parameters. 
The  even  moments  of  the  weights  (summed  over  the  fitting  times  {/.,;})  would  then  be  the  [n  -  1] 
non-linear  search  parameters.  If  the  problem  is  intrinsically  asymmetric,  then  there  would  be 
[2n  —  2]  moments  to  use  as  nonlinear  search  parameters.  For  extrapolation,  it  seems  clear  that 
there  will  be  little  likelihood  of  driving  any  negative,  but  it  remains  a  ccmcern  for  the  general 
case  and  must  be  guarded  against. 

(Consider  for  example  the  case  of  choosing  a  weighted  least-scjuares  fit  of  a  general  quadratic  to  N 
phase  comparison  data  points  at  a  specific  set  of  times  {f,}.  For  a  specific  noise  model  described 
by  the  coefficients  {/ia})  we  want  to  choose  the  weights  to  minimize  the  standard  uncertainty  in 
the  average  frequency  over  the  time  interval  [f ,  t  -I-  t]  .  By  constraining  weight.s  to  sum  to  1 ,  and 
by  constraining  the  first  through  fourth  moments  of  the  weights  to  be  independent  of  the  first 
[N  “  4]  weights,  we  can  ensure  that  G~^  is  independent  of  [N  —  4]  w'eights.  By  equating  to  zero  the 
[N  —  4]  partial  derivatives  of  u.^(Tt)  with  respect  to  VF,  we  can  minimize  the  standard  uncertainty 
in  average  frequency  with  respect  to  these  [jV  —  4]  weights.  The  easiest  form  to  differentiate  for 
this  purpose  is  one  like  that  of  Equation  8,  which  has  collected  all  the  terms  multiplied  by  any 
weight  Wi-  Including  the  constraint  equations,  we  then  have  N  linear  equations  in  the  N  unknown 
weights  {kFi},  parameterized  in  the  4  moments  remaining  to  be  searched.  The  optimized  standard 
uncertainty  for  this  set  of  four  moments  is  evaluated,  and  a  four- parameter  search  (ea('h  set  of 
moments  being  optimized  by  re-solving  the  N  linear  equatiorcs)  this  search  is  tractable  by  the 
simplex  method  (for  example).  If  the  problenv  is  symmetric  about  some  time  (symmetry  for  both 
{/j.}  and  [/-,f  -F  rj),  it  can  be  set  up  so  that  the  first  and  third  moments  are  zero,  and  there  would 
l>e  only  two  parameters  to  search, 

Choasing  weights  is  simpler  for  a  linear  least-squares  fit  to  N  pha.se  compari.son  data  points, 
taken  at  a  specific  set  of  times  {t?,}.  To  metafit  the  best  weights  that  minimize  the  standard 
uncertainty  in  the  average  frequency  over  the  interval  [t,t  +  r]  for  the  noise  model  of  interest, 
described  by  the  coefficients  {ha}i  we  can  again  linearize  the  problem  -  but  with  only  two  search 
parameters  (the  first  and  second  moments  of  the  weights).  We  define  three  constraint  equations 
Sfli  =  Ail  =  M2.  The  N  partial  derivatives,  with  respect  to 

the  weights,  of  the  standard  uncertainty  in  average  frequency  over  the  interval  [t,  t,  -F  r]  give  a  set 
of  N  equations  F  •  IF  =  f,  where  Fij  =  {t,,  —  Mi)t/{M2  —  Mf)  <  [x()(ti)  —  a;o(fi)][a;o(f,)  —  2:o(fi)]  > 
and  r  j  =<  [xQ{t  -Ft)  —  a;o(f)][2Jo(f ,)  —  2:o(fi)]  >.  The  finst  column  of  F  is  a  column  of  zeros.  Three 
of  these  equations  are  to  be  replaced  by  the  three  constraint  equations;  one  replacement  is  for  the 
most  ill-conditioned  equation  j  which  has  tj  clo.sest  to  the  centroid  of  the  weights  (Mi)  for  this 
iteration,  the  other  two  replacements  2ire  more  arbitrary.  If  the  problem  is  symmetric  about  some 
time  (symmetry  for  both  {f,}  and  [f,  f  -F  r]),  it  can  be  set  up  so  that  the  first  moment  is  zero,  and 
there  would  be  only  one  parameter  to  sear  ch. 
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An  even  simpler  case  of  metafil.ting  is  the  choice  of  weights  in  a  simple  w’eighted  average,  for  multiple 
calibration  runs  to  minimize  the  standard  uncertainty  in  the  average  frequency  for  a  specific  period, 
arbitrarily  placed  with  respect  to  the  calibration  runs.  We  consider  calibration  intervals  long  enough 
to  be  in  the  regime  where  the  tw’o  end  point  method  is  chosen  for  each  calibration  run,  with  M 
such  calibration  intervals  [tj,  ,  + Tj.].  For  the  weighted  average  of  the  M  calibrations,  the  standard 
uncertainty  in  the  average  frequency  over  an  interval  [t,t.  -1-  r],u^(t,T)  is 


Assigning  a  weight  of  —1  to  the  interval  [t,  t  +  r],  defining  tq  as  being  equal  to  r,  Equation  11  can 
be  rewritten  as 


t  +  r) 


(12) 


A  solution  for  the  optimum  weighting  procedure  is  relatively  easy  to  find  since  the  minimum  value 
for  Uy{t,  t.  +  t)  is  to  be  found  for  values  of  tn,  satisfying  t  +  t)]  =  0,  so  that  after  taking 

the  derivative  and  separating  out  the  i  =  0  term 


^  +  T,:)  -  2;(t|)][x(tfe  +  Tfc)  -  x{tk)]) 

^  — ([ic(t  +  r)  -  x{t)][x{tk  +  Tk)  -  s(tfc)]).  (13) 

T'k 

We  use  M  —  1  of  these  equations,  and  for  the  equation  we  use  the  normalization  equation 
of  the  weights:  Wj,  —  1.  This  gives  M  simultaneous  linear  equations  in  the  M  unknown 

weights.  The  general  interval  covariance  has  analytic  forms  for  our  noise  model,  in  terms  of  the 
Z-function  [4],  If  we  define  the  M  x  M  matrix  F:  Fj  ,  =  1  for  j  =  I..M,  -f  r,  —  tj) 

+  2{tj  +Tj  —  tj)  -  Z(b;  +  Ti  —  tj  ~  Tj)  -  Z{t,:  —  tj)]  for  ?.  =  2..M  and  j  =  1..A/,  and  define  f:  t'l  =  1 
and  Tj  =  /?‘aclrTj[J(t  +  t  —  tj)  +  2(tj  +  Tj  —  t)  -  Z(t  +  r  —  tj  —  Tj)  -  I(t  -  t,)]  for  j  =  2..M.  The 
M  dimensional  weights  vector  riJ  is  F~F 


Applications 


For  any  given  potential  application  of  metafitting  weights,  we  must  con.sider  whether  metafitting 
is  more  than  an  interesting  academic  exercise;  can  metafitting  find  a  reduction  in  the  standard 
uncertainty  which  is  a  significant  improvement?  Since  uncertainties  are  rarely  established  to  better 
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than  10%,  an  improvement  should  be  larger  than  this  to  be  deemed  significant.  Therefore  we  have 
examined  the  simplest  case,  of  linear  extrapolation,  discussed  above,  and  for  the  five  different  power- 
law  noise  types  we  have  considered  distributions  of  weights  with  different  moments  [6],  We  have 
examined  the  expected  standard  uncertainty  for  both  symmetric  extrapolation  suited  to  time-.sca]e 
calibration  (where  post- processing  can  be  used  to  apply  calibrations  from  the  “future”)  and  to  time- 
assymetric  extrapolation  suited  to  real-time  applications.  For  .symmetric  extrapolation  intervals 
that  are  large  compared  to  the  calibration  run’s  duration,  different  common  weight  distributions 
gave  similar  uncertainties  (differing  by  less  than  10%)  except  for  white  phase  noise.  For  one-way 
extrapolation  for  time.s  much  longer  than  the  calibration  run’s  duration,  the  uncertainties  are  even 
more  similar  (less  than  2%  advantage  for  end-point  fitting  over  equal  weights,  except  for  white 
phase  and  flicker  phase  noise).  Thus  for  many  PTTl  applications,  end-point  fitting  and  equal- 
weight  fitting  give  similar  standard  uncertainties,  and  the  choice  should  be  between  the  greater 
simplicity  of  the  end-point  fit  and  the  greater  robustness  of  the  equal-weights  fitting  procedure. 

In  real-life  PTTI  work,  robustnass  would  often  prevail  over  simplicity.  For  trying  to  optimize 
results  from  multiple  calibration  runs,  simplicity  is  valuable  to  us  while  robustness  is  not  needed 
in  the  model.  The  optimum  processing  of  a  number  of  calibration  runs  is  expected  to  be  largely 
independent  of  the  processing  within  the  run. 

The  main  application  which  has  attracted  our  attention  is  the  optimal  use  of  hydrogen  masers, 
calibrated  periodically  in  frequency  with  intermittently  operated  cesium  fountain  frequency  stan¬ 
dards  [8],  [6].  We  consider  two  type.s  of  maser  operation:  free-running  and  autotuned.  We  use  two 
power  law  models  for  the  maser  noise,  representing  a  free-running  hydrogen  maser  (type  1)  with 
/i2  =  2.7x10-24^ /ij  =  2.9x10-30,  ho  =  2.9xl0-2^fi_i  =2.6x10-31  and /i_2  =  7.2 x  10-3®;  and  an 
auto-tuned  maser  (type  2)  with  I12  =  6.7xl0-^3^  =  2.9x10-30,  ho  =  2.9x10”^^,  =  7.2x10-31 

and  fi_2  =  4.9  x  10-3'’.  has  two  low-flux  masers  which  would  benefit  from  a  metafitting 

optimization  of  the  weights  within  a  calibration  run  of  an  hour,  since  there  is  still  some  white  phase 
noise  contribution  for  this  calil)ration  interval.  Preliminary  analysis  suggasts  that  the  end-point 
procedure  is  within  10%  of  the  optimum.  For  phase  data  taken  every  30  s  for  an  hour,  extrapolated 
to  an  interval  of  a  day,  the  end  point  method  is  1.2%  better  than  the  equal-weight  linear  least 
squares  fit  for  our  free- running  maser  model,  and  as  good  for  the  type  2  maser  model.  Thus  we 
can  use  the  simple  two  end  points  procedure  to  establish  the  l?est  frequency  transfer  accuracy  for 
multiple  calibration  runs.  For  this  procedure  the  standard  uncertainty  for  multiple  calibration  runs 
can  be  calculated  more  easily  than  in  the  general  case. 

Within  the  context  of  end-point  fitting  from  each  calibration  run  there  are  still  metafitting  choices 
to  be  made  about  the  way  in  which  the  runs  are  to  be  used.  One  po.s.sib]e  strategy  is  a  loose  lock  in 
frequency:  after  a  calibration  run  (an  hour  in  duration,  in  our  example)  is  complete,  the  frequency 
of  the  maser  is  reset  (through  the  synthesizer  control,  for  example),  either  immediately  -  or  after 
some  delay.  Clearly  the  least  delay  is  best,  and  we  chose  this  procedure  with  zero  delay  as  the 
reference  procedure  as  we  examine  a  series  of  possilfie  improvements. 

A  slightly  better  possibility  might  be  to  have  an  output  tightly  locked  in  phase  to  the  cesium 
fountain  during  the  calibration  run,  followed  by  a  frequency  lock  to  the  fitted  frequency  of  the 
calibration  run.  The  phase-lock  type  of  frequency  control  removes  the  noise  of  the  maser  during 
the  calibration  run,  giving  it  an  advantage  that  remains  noticeable  for  extrapolation  intervals  many 
times  longer  than  the  calibration  interval.  However,  for  extrapolations  of  an  hour-long  calibration 
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out  to  a  period  of  a  day  or  more,  there  is  not  a  large  advantage:  2.3%  for  the  free-runniiig  maser 
and  2.4%  for  the  autotuned  maser  model. 

A  more  .significant  advantage  comes  from  allowing  postprocessing,  as  can  often  be  tolerated  in  time- 
scale  construction  and  for  frequency  intercomparisons.  We  consider  a  single  calibration  interv'al 
t-c  and  calculate  the  ratio  of  the  standard  uncertainty  of  the  average  frequency  over  an  interval  r 
for  the  best  real-time  frequency  control  to  the  .symmetrically  extrapolated  time  interval  r.  The 
quantitative  postprocessing  advantage  will  depend  upon  the  specific  processing  scheme  or  schemes 
envisaged  -  the  duration  and  frequency  of  calibration  intervals.  The  po.stprocessing  advantage  is 
up  to  a  factor  of  two  [6]. 

A  postproceasing  advantage  of  two  is  really  quite  significant.  To  achieve  the  same  improvement 
in  the  maser  ensemble  could  be  done  -  by  increasing  the  maser  ensemble  size  by  four  times.  The 
postprocessing  advantage  of  greatest  interest  to  us  is  for  r  repre.senting  extrapolation  to  the  time 
interval  between  calibrations  -  which  we  expect  would  l)e  between  1  day  and  1  week.  Initial  interlal> 
oratory  frequency  intercomparisons  between  cesium  fountains,  before  regular  calibration  schedules 
can  be  set  up,  may  require  extrapolation  times  longer  than  1  week  for  minimum  uncertainty. 

Envisaging  multiple  frequency  calibration  runs  per  week,  of  either  hydrogen  maser  type  with  a 
cesium  fountain  having  a  standard  uncertainty  of  optimistically  5  per  week,  at  the  same 

time  each  working  day,  what  is  the  best  weighting  procedure  for  using  these  calibrations  in  an 
algorithm  to  determine  the  frequency  over  a  given  interval?  For  the  week’s  pattern,  postprocessing 
extrapolation  of  each  day’s  results  independently,  using  the  frequency  from  the  nearest  calibration 
interval  gives  a  77%  improvement  in  accuracy  for  the  free-running  maser,  and  an  improvement  of 
29%  for  the  auto- tuned  maser. 

We  have  solved  for  the  optimum  weights  of  the  maser  calibrations  to  give  the  lowest  standard 
uncertainty  in  average  frequency  over  one  week  [6].  The  week  is  best  spanned  by  weighting  Monday 
and  Friday  runs  more  heavily,  to  account  for  the  weekend  gap  in  calibration.s.  For  the  type  1  maser, 
the  optimum  weights  follow  the  spanning  times  rather  closely,  and  the  optimum  weights  offer  only 
a  1.1%  improvement  in  average  frequency.  For  a  type  2  maser,  there  is  a  4.7%  improvement. 

If  adj accent  weeks’  calibration  runs  are  also  available,  and  the  average  frequency  over  a  particular 
week  is  required,  the  optimum  metafitting  includes  a  small  admixture  from  the  preceding  and  the 
following  weeks.  For  a  type  1  maser,  most  of  the  weight  comes  from  the  preceeding  Friday  and  the 
following  Monday.  For  an  autotuned  (type  2)  hydrogen  maser  noise  model,  the  optimum  weights 
have  a  .slower  variation  through  the  weeks,  and  the  three-week  optimum  has  several  %  of  the  weight 
on  point.s  that  are  a  full  week  from  the  calibration  rums  of  the  central  week.  There  is  a  19% 
improvement  to  the  t3'pe  1  maser,  and  an  18%  improvement  for  the  type  2  maser.  The  improvements 
are  summarized  in  Table  I,  given  with  standard  uncertainties  and  cumulative  advantages  as  each 
improvement  is  applied.  For  either  maser  model,  the  optimization  of  weights  to  apply  to  each  run 
over  multiple  weeks  gives  about  a  20%  improvement  in  accuracy  from  the  equal-weight  ca.se.  It  i.s 
not  a  large  improvement,  but  it  is  almost  free  -  although  it  does  give  additional  cross-correlation 
between  each  week’s  frequency  processed  in  thi.s  way.  Ca.scaded  with  the  other  advantages  discussed 
earlier,  it  results  in  a  factor  of  2.2  improvement  in  the  accuracy  transferrable  with  a  free-running 
(type  1)  hydrogen  maser;  and  an  improvement  of  64%  for  the  auto-tuned  (type  2)  maser.  For  the 
free-running  maser  model,  the  inetafitted  optimurrt  standard  uncertainty  is  6.8  times  sirtalier  than 
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Type  I 

Cum. 

Type  2 

Cum. 

method 

Uy{7d) 

Adv, 

adv. 

Adv. 

adv. 

I 

/  reset  to  unw'eighted  fit 

1.70  X  10-'''’ 

DSSSSIfl 

mwm 

II 

/  reset  to  end  points 

1.77  X  lO-'-'' 

1.2% 

■KQIl 

III 

phase  lock  -f  II 

T73  X  10-''^ 

2.3% 

1.04 

TT2  X  10-*^ 

2.4% 

1.02 

IV 

daily  postprocessed 

I’flilffftViiM 

77% 

1.83 

29% 

V 

metafit  1  week 

1.85 

4.7% 

1.39 

VI 

metafit  3  weeks 

19% 

2.21 

msm 

VII 

metafit  5  weeks 
.. 

0% 

1.64 

CTy{7d) 

5,48  X 

T72  X  10-'''' 

_ _ _ 

Table  1 :  Reduction  of  standard  uncertainty  in  average  frequency  at  7  days,  for  a  free-running  (type 
1)  maser,  and  an  autotuned  maser  (type  2),  when  controlled  by  different  methods  from  five  1-hour 
calibrations  per  week.  The  %  advantage  for  each  method  is  the  accuracy  improvement  over  the 
previous  method.  The  last  column  gives  each  method’s  cumulative  advantage  over  method  I,  a 
synthesizer  reset  to  the  least-squares  calibration  fit.  The  Allan  deviation  (Ty(T  =  7d)  is  also  given. 


the  Allan  deviation  at  1  week,  and  for  the  type  2  maser  it  is  2.5  times  smaller  than  the  Allan 
deviation  at  1  week. 

Other  interesting  strategies  are  beyond  the  scope  of  this  w'ork.  Longer  runs  on  Monday  and  Friday 
and/or  early- Monday  and  late-Friday  calibration  runs  could  be  invoked  to  further  improve  the 
performance.  Our  methods  allow  for  weight  optimization  for  any  set  of  calibration  runs,  and  for 
calculating  the  resulting  standard  uncertainty  in  average  frequency, 

For  some  applications,  statistical  independence  of  each  week,  or  each  10-day  period,  may  be  highly 
valued  -  for  example,  the  clock  reports  to  BIPM  each  10  days  that  are  used  for  determining  TAI 
(and  UTC)  should  be  independent  of  each  other.  Weights  for  data  from  the  weekly  caliijration 
cycle  could  be  re-optimized  for  the  seven  different  10-day  cycles  that  would  exist.  The  metafitted 
optimum  weights  for  the  two  maser  models  are  shown  in  Figure  1.  For  the  free-running  maser 
model,  the  70-day  standard  uncertainty  in  average  frequency  is  3.31  x  10“^®  for  the  combination  of 
the  seven  independent  optimized  10-day  periods,  as  compared  to  3.07  x  10“^**  for  the  combination 
of  10  independent  7-day  periods.  For  the  autotuning  maser  model,  the  70  day  standard  uncertainty 
in  average  frequency  is  2.72  x  10“^®  for  the  combination  of  the  seven  independent  optimized  10-day 
periods,  as  compared  to  2.62  x  10~^°  for  the  combination  of  10  independent  7-day  periods. 


Conclusion 

Our  method  for  calculating  the  standard  uncertainty  for  realistic  noise  models  has  allowed  us  to 
compare  a  wide  variety  of  algorithms  for  treating  one  particular  calibration  schedule.  We  have 
metafitted  the  algorithm  in  several  ways,  and  have  identified  ways  to  improve  the  accuracy  of  the 
maser  frequency  control  by  2.2  and  1.64  times.  We  find  that  using  the  10-day  BIPM  schedule,  with 
independent  processing  of  the  calibrations  for  the  10-day  periods,  the  expected  asymptote  for  a 


358 


single  auto-timed  (type  2)  maser  could  reach  1 .2  x  10“^'^  at  1  year.  For  a  fre(^niniiing  (type  1)  maser, 
the  standard  uncertainty  at  1  year  would  be  1.5  x  10“^^  .  Thus  a  flicker  floor  and  accuracy  of  10“^® 
for  the  cesium  fountain  is  accessible  for  periods  of  a  year  with  current  masers  carrying  the  time 
scale.  Operating  the  masers  at  the  stability  level  of  the  masers  presents  a  challenge.  Transferring 
10“^^’  frequency  accuracy  to  a  second  laboratory  also  pre.sents  a  challenge.  The  reliability  of  a 
cesium  fountain  which  might  do  this  seems  to  be  a  major  challenge,  perhaps  comparable  to  the 
challenge  of  making  a  cesium  fountain  with  a  flicker  floor  and  accuracy  of  10“^^.  Perhaps  the 
greatest  value  of  thi.s  metafitting  procedure  i.s  to  .show  the  very  best  performance  which  might  be 
extracted  from  masers  represented  by  these  models.  If  greater  accuracy  is  desired,  then  different 
approaches  must  be  used. 
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Autotuned  Maser 
Free-running  Maser 


Optimum  Weights  for  Weekday  Calibration 
of  10-day  Average  Frequency 


Figure  1.  Optimum  weights  for  combining  weekday  calibrations  that  give  the  minimum  standard 
uncertainty  in  the  average  frequency  over  a  10-day  interval.  The  calibrating  reference  standard  is  taken  to  be  an 
ideal  one,  used  for  one  hour,  at  the  same  time  every  working  day.  The  optimum  weights  are  shown  for  two  flywheel 
oscillators:  a  free-running  hydrogen  maser  model  and  an  auto-tuned  maser  model.  The  optimum  weights  are  shown 
for  a  the  10-day  period  starting  on  each  day  of  the  week. 
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